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Abstract. We present a simple method for the solution of one-compo- 
nent and multicomponent hydrodynamic equations based on the Newton- 
Raphson method. We show that this method can be used for the solution 
of stationary hydrodynamic equations. This method has been used for 
the calculation of the low density stellar wind models for which the mul- 
ticomponent nature of the wind influences the overall wind structure. 



1. Introduction 

Radiatively driven stellar winds are accelerated by the absorption of radiation 
mainly in the resonance lines of such elements like C, N, O or Fe (see Castor, Ab- 
bott, & Klein 1975, hereafter CAK; Abbott 1982). These wind components have 
much lower density than the stellar wind itself, which is composed mainly by hy- 
drogen and helium. The momentum obtained by the low-density absorbing wind 
component is transfered to the high-density nonabsorbing wind component via 
collisions between electrically charged particles (Springmann & Pauldrach 1992; 
Babel 1995). Clearly, multicomponent radiatively driven winds of hot stars have 
a multicomponent nature. This multicomponent nature is important mainly for 
the low-density stellar winds, which shall be described by the multicomponent 
hydrodynamic equations. 

Krticka & Kubat (2000; 2001a; 2001b) calculated models of such multicom- 
ponent stellar winds using Newton-Raphson method. In this paper we describe 
in detail method of solution of these multicomponent hydrodynamic equations. 



2. Model equations 

Each wind component is described by a density p a , radial velocity v ra , tempera- 
ture T a and electrical charge q a . Typically, we take into account only three wind 
components (ie. absorbing ions, nonabsorbing ions and electrons). 

Stationary spherically symmetric multicomponent radiatively driven stellar 
winds are described by the continuity equation for each component a of the wind 

(r 2 p a Vra) = 0, (1) 

(note, that in the case of continuity equation of electrons a right-hand term 
occurs due to the change of the wind ionization in the flow, see Krticka & 
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Kubat, 2001b), momentum and energy equations for each component a 



where a a is isothermal sound speed, E is charge separation electric field, g is 
gravitational acceleration, radiative acceleration g™ d acts on free electrons and 
absorbing ions. Radiative acceleration acting on absorbing ions is taken in the 
CAK approximation with the finite disk correction factor (Pauldrach, Puis, & 
Kudritzki 1986) and with force multipliers after Abbott (1982), slightly modified 
for the multicomponent case (Krticka & Kubat 2000). Frictional acceleration 
g^°, frictional heating Q^ c and exchange of the heat between wind compo- 
nents a and b are given by Burgers (1969; see also Krticka & Kubat 2001b). 

The only transitions which contribute to the radiative heating/cooling term 
Qrad - m static atmosphere are bound-free and free-free transitions (Kubat, 
Puis, & Pauldrach 1999; see also Kubat, this volume). Because these transitions 
deposit energy directly to electrons, this radiative heating/cooling term is con- 
sidered in electron energy equation. Mean intensity J u at the base of the wind 
is taken as an emergent radiation from a spherically symmetric static hydrogen 
model atmosphere (Kubat 2001). In the case of moving media additional heat- 
ing/cooling term occurs, so-called Gayley-Owocki heating, which is caused by 
the dependence of a radiative force on a velocity via Doppler effect (Gayley & 
Owocki 1994; Krticka & Kubat 2001b). 

The electrical charge q a is obtained using the approximate ionization equi- 
librium described by Mihalas (1978, eq. [5-46]). 

The system of hydrodynamic equations (1, 2, 3) is closed using an equation 
for polarisation electric field in the form of 



Finally, the system of equations solved is supplemented by appropriate 
boundary conditions. We start to calculate our models at the critical point 
of nonabsorbing component. We also use zero current condition, the condition 
of quasi-neutrality and assume that the flow at the inner boundary is in radiative 
equilibrium and that the boundary temperature of all components is the same. 
Boundary density of absorbing ions is determined using abundance ratios. 

3. Method of solution 




(2) 



(3) 





Velocities, densities and other unknown variables are calculated at grid points. 
We selected logarithmically spaced grid in radii with accumulation of grid points 
near the stellar photosphere, 
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where i is a number of a given grid point and NR is number of grid points. We 
use typically 50-300 grid points. 

For the approximation of the derivative of expression X we selected 



dX 
dr 



i = NR, 



(6) 



where 



Ar i = r i -r i - 1 , n = - (r { + rj_i) , Vi = zr 1 — (7) 

2 r i+ i - rj 

An appripriate approximation of the derivative is important, numerical tests 
showed that approximation (6) gives the best convergence of the models. 

An approximation of hydro dynamical equations (1, 2, 3) together with the 
equation for a polarisation electric field (4), equation of ionization equilibrium, 
and boundary conditions can be formally written as 

P</> = 0, (8) 
where the vector describing the solution has the form of 

rj) = (V>i,i/>2>--- >V>nr) T > (9a) 

\ a a a a / 

and P consists from all equations mentioned above. 

For the solution of the nonlinear system of model equations (8) we selected 
Newton- Raphson method. The solution of Eqs.(8) can be obtained using itera- 
tive procedure in the form of 

rSrp n+1 = -P> n , (10) 

where tj) n denotes solution in the n-th iteration, 8tl) n+l is a correction of the 
solution and the Jacobian is 

J m ~ q^- (ii) 

The Newton- Raphson method is similar to the well-known method of complete 
linearization (Auer &; Mihalas 1969). The analytical form of Jacobi matrix 

can be easily obtained from Eqs.(8) and for the case of three-component 
radiatively driven stellar wind was calculated by Krticka (2001). However, the 
value of corrections Si^ n+1 in each iterative step is limited by 0.7ip n during 
the calculation of a model below the critical point and by 0.lip n during the 
calculation of a model above the critical point. The velocity differences between 
wind components are limited similarly. 

Another problem for the calculation of wind models is the inclusion of crit- 
ical point condition. Critical point is a point where the wind velocity is equal to 
the sound speed. Generally, there are two possibilities of inclusion of the critical 
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point condition. Nobili & Turola (1988) showed that critical point condition 
can be directly included into the set of linearised equations. However, due to 
some numerical problems, we selected the so called "shooting method". We 
change the density at the base of the wind in order to fulfil the critical point 
condition. When the base density and model upstream the critical point are 
obtained, we calculate model downstream the critical point. The initial model 
for the Newton-Raphson iterations is selected in such a way that it converges 
to the correct branch of CAK solution appropriate for the flow downstream the 
critical point (see CAK). 

For the solution of the system of linearised equations (10) we use numerical 
package LAPACK. We also applied the Gaussian elimination, however we found 
that the package LAPACK is faster and even more accurate. 

Newton-Raphson method represents quite powerful tool for the solution of 
stationary hydrodynamic equations. The iterations converge very fast (typical 
relative change is of order 10~ 10 after 20 iterations), the number of equations 
solved can be easily extended and it is possible to calculate multidimensional 
models. On the other hand the model convergence depends on initial estimate 
and we should also check the dynamical stability of obtained solution. 
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